tos o. -os- iTSr.cf^ 


NASA Contractor Report 178028 

ICASE REPORT NO. 85-60 NASA-CR- 178028 

19860009180 

ICASE 

MARCHING ITERATIVE METHODS FOR THE PARABOLIZED 
AND THIN LAYER NAVIER-STOKES EQUATIONS 


Moshe Israeli 


Contract No. NAS1-17070 
December 1985 


INSTITUTE FOR COMPUTER APPLICATIONS IN SCIENCE AND ENGINEERING 
NASA Langley Research Center, Hampton, Virginia 23665 


Operated by the Universities Space Research Association 


NASA 

National Aeronautics and 
Space Administration 

Langley Research Center 

Hampton. Virginia 23665 



3 1176 01308 7201 

ns 



MARCHING ITERATIVE METHODS 

FOR THE PARABOLIZED AND THIN LAYER NAV IER-STOKE S EQUATIONS 


Moshe Israeli 

Technion - Israel Institute of Technology, Haifa (Israel) 


Abstract 

Downstream marching iterative schemes for the solution of the 
Parabolized or Thin Layer (PNS or TL) Navier-Stokes equations are 
described. Modifications of the primitive equation global relaxation 
sweep procedure result in efficient second-order marching schemes. 

These schemes take full account of the reduced order of the approximate 
equations as they behave like the SLOR for a single elliptic equation. 
The improved smoothing properties permit the introduction of Multi-Grid 
acceleration. The proposed algorithm is essentially Reynolds number 
independent and therefore can be applied to the solution of the subsonic 
Euler equations. The convergence rates are similar to those obtained by 
the Multi-Grid solution of a single elliptic equation; the storage is 
also comparable as only the pressure has to be stored on all levels. 
Extensions to three-dimensional and compressible subsonic flows are 
discussed. Numerical results are presented. 
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1* INTRODUCTION 


Considerable evidence accumulated recently about the applicability 
of the Parabolized Navier-Stokes equations for high Reynolds number flows 
with a principal flow direction; see Rubin [1] . The PNS equations are 
obtained by neglecting the streamwise viscous terms in the Navier-Stokes 
(NS) equations. When the viscous terms in the circumferential direction 
are also neglected, one gets the Thin Layer approximation. 

The steady PNS equations still have an elliptic nature, and there- 
fore the initial value problem in the marching direction is not well 
posed [2] . A well posed initial-boundary value problem can be formulat- 
ed by specifying (for example) upstream and side conditions for the 
velocities and one downstream condition for the pressure. This coupled 
system of partial differential equations behaves like a single elliptic 
equation for the pressure. Therefore the PNS equations must be solved 
globally and cannot be solved by a single sweep marching. The reduced 
order of the PNS equation can be exploited by constructing an iterative 
marching method for updating the pressure field only. Such a multiple 
sweep iteration method has the advantage that the velocity fields are 
generated during the marching process and only the pressure field has to 
be stored from sweep to sweep. A considerable saving in storage results. 
However, simple minded marching does not result in good convergence 
properties and sometimes diverges. For the two-dimensional incompres- 
sible case, Israeli and Lin [3] devised a stable marching scheme that 
behaves like the Successive Line Over Relaxation (SLOR) method for a 
single elliptic equation. The good smoothing properties of the above 
mentioned scheme can be used in a Multi-Grid (MG) framework in order to 
accelerate the convergence of the solution of the PNS (or TL) equations. 
The marching scheme is implemented using a new stable algorithm which is 
second order also in the marching direction. The same method can be 
used without modification for the subsonic Euler equations as the effect 
of the Reynolds number on the convergence rate is insignificant. In two 
dimensions the PNS and TL equations are identical and therefore the same 
analysis applies to both. 
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It turns out that the extension to three-dimensions is conceptually 

simple; but the resulting algorithm, a successive plane over relaxation, 
is complicated by the requirement of the simultaneous solution of the 
equations in planes perpendicular to the marching direction. This 
problem can be alleviated by splitting of the equation of continuity 
from the momentum equations. 

The extension of the method to compressible flows is conceptually 
non trivial. The original iterative method is based on the concept 
that the convergence relies on the implicit relaxation of a single 
quantity, the pressure, which approximately satisfies a single elliptic 
equation. In the compressible case a viable approach is to eliminate 
the pressure and to derive an equation for p, the logarithm of the 
density. It can be shown that p satisfies approximately equation (1.1) , 


(1 - + 


ss 


P = o, 
*nn 


( 1 . 1 ) 


where M is the Mach number and s and n are coordinates along and 
perpendicular to the flow direction. Although this equation is never 
derived or used in the algorithm, it reveals the fact that for M < 1 
the upstream influence is transmitted through the quantity p, and 
therefore only this quantity should be stored or updated. The flow of 
information should be downstream for the velocity and temperature and 
upstream for the density, and the difference scheme must be built 
accordingly. For supersonic flows the flow of information should be 
only downstream and the marching method is non-iterative. For super- 
sonic flows with imbedded subsonic regions, the iterative method should 
be used, combined with an appropriate switching at shock waves and sonic 
lines . 

It should be pointed out here that the present approach is very 
different conceptually from that of Reddy and Rubin [4]. Although they 
used our idea* of backshifting the pressure, one full mesh distance, 
with respect to the velocity for incompressible flows, their 
generalization to compressible flows is a Mach number dependent shift 


^Israeli, M. (1982), NASA Lewis seminar, July 1982. 
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which vanishes for M > 1. This smooth transition from subsonic to 
supersonic flows is questionable since the change of type of equation 
( 1 * 1 ) is sudden, at M = 1 . Indeed, only our full shift is used in their 
papers and properly results in a conservative scheme across a shock. 

Another question raised by the above mentioned paper is that of the 
distinction between the pressure which uses downstream data and the 
density which uses upstream data. This obscures the issue of the direc- 
tion of flow of information and proper location of boundary conditions. 
This approach should result in inconsistency of boundary data and may 
eventually lead to ill posedness and divergence. 

In the next sections we will summarize our previous theoretical 
results, present some new numerical results and the extensions to 3-D 
and compressible flows. 


2. FORMULATION FOR THE INCOMPRESSIBLE CASE 

For simplicity we will consider initially the case of the steady, 
incompressible, and two-dimensional PNS (or TL) equations in cartesian 
coordinates [x;y] : 

U + V = 0 
x y 

(u 2 ) + (UV) = -P + U /Re 
x y x yy 

(UV) + (V 2 ) = -P + V /Re 

x y y yy 

where x is the mainstream direction. Re is the Reynolds number, 
and V are the nondimensional velocity components in the x and 
direction, respectively. P is the nondimensional pressure. 

The two-dimensional NS equations are elliptic of order four - Brandt 
and Dinar 1 5 J . The PNS are elliptic only of order two like the Poisson 
equation (the mathematical nature of several two-dimensional and three- 
dimensional approximations to the Navier-Stokes equations was analyzed 
in [6]). This ellipticity is due to the pressure gradient terms via 
the continuity equation. A well posed problem can be formulated by 
defining the boundary conditions as described in Fig. 2. The following 
Dirichlet conditions may be specified: 


( 2 . 1 ) 

( 2 . 2 ) 

(2.3) 

U 

y 
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upstream boundary (AB) : 

U = U. 

m 

V = V. 
in 

(2.4) 

at a solid wall (AD) : 

U = U 

wall 

V = V 

wall 

(2.5) 

at the outer boundary (BC) : 

U = U . ? 

out 

V = V 

out 

(2.6) 

at the downstream boundary (CD) : 

P = P , 

down 


(2.7) 


Other boundary conditions can be used, but the same number of conditions 
on each boundary must be kept. 

In order to separate linear and non linear effects, some of the 
convergence tests were performed with the following linear version of 
equations (2,1)~(2.3): 


U + V = 0 
x y 


(aU) + (bU) = -P + U /Re 

x y x yy 

( aV ) + (bV) = -P + V /Re 

x y y yy 


where a and b are known functions of x and y. 

3. DISCRETIZATION AND MARCHING 

Numerical solutions of Eqs. (2.1)- (2.3) are obtained by spreading 
a grid over the computational domain. Let us assume that the grid points 
are distributed evenly along the x and y coordinates with the spac- 
ing Ax and Ay respectively. When differencing these equations it 
should be remembered that their nature should be reflected [1,7] in 
the finite difference approximation. In order to be consistent with the 
boundary layer (parabolic) nature of the flow, the axial gradients of 
the velocities should be computed using only upstream values, while the 
elliptic nature is preserved by forward differencing the axial pressure 
gradient [1,7,8]* Consequently, it was assumed that a stable marching 

scheme must be of the first order in the marching direction. It turns 
out that this effect can be achieved by a judicious choice of the place- 
ment of the variables to be solved at each station. The choice can be 

explained most easily by taking V = 0 and — = 0 in Eq. (2.2) for U, 

Re 

yielding 


( 2 . 8 ) 

(2.9) 

( 2 . 10 ) 
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A first order difference scheme then becomes 
2 2 

m, 3 m-1,3 m, 3 ^m+1,3 * 

the unknowns are U . and p , . The scheme first suqqested by Israeli 
mo m, 3 ” * 

[ 9 , 10 ] is: 

2 2 

U , . - U . = p . - p , 

m+l,j m, 3 m, 3 m+l,j 

with the unknowns U _ . and p . . The scheme is centered about m + — 

m+1,3 m, 3 2 

and is second order. This approach was subsequently used by Rubin and 
Reddy [8] and Reddy and Rubin [ 4 ] . 

In addition, one may stagger the velocity V with respect to the 
other variables as shown in Fig. 3 , where the centering points of the 
different difference equations are also plotted. The differential equa- 
tions are approximated by central second-order approximations whenever 
needed averaging was used as is usually done for staggered grids. 

Numerical experiments with a first order computer code show that 
the solution after one marching sweep is not close to the final solution 
of the PNS equations when the initial pressure field is constructed 
using the boundary layer assumption p^ = 0 . Since the p^ term is 
forward differenced, some global iterations over the whole solution 
domain should be performed in order to converge the explicit contribu- 
tion to this pressure terp. The simplest global iterative technique to 
solve the equations is by multiple marching sweeps with the primitive 
equations where only the pressure field is kept from iteration to itera- 
tion [ 1 ] . Numerical experiments also show that for certain nets this 
procedure diverges. The divergence occurs also for the linearized ver- 
sion of Eqs. ( 2 . 1 ) - ( 2 . 3 ). Figure 1 presents the residual of the pressure 
field as a function of the global iteration's sweep number for a 21 x 11 
field. A jump is encountered every 10 iterations (probably related to 
the arrival of the boundary pressure pulse traveling at the numerical 
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scheme speed) leading to ultimate divergence. However convergence was 
reported with different mesh and boundary conditions and also when 
combining the above procedure with a multigrid technique [4] . It was 
thought that the replacement of one of the momentum equations by the 
Poisson equation for the pressure will improve the convergence rate, but 
the solution did not satisfy the replaced momentum equation. A success- 
ful implementation of the marching technique is derived in the next 
section. A short and reduced version of the analysis was presented first 
in [3]. 


4. A MULTI-GRID ALGORITHM 

The Multi-Grid technique is a numerical strategy for substantially 

improving the convergence rate of an iterative procedure. In order to 

facilitate comparison with theory , the accomodative C-cycle MG algorithm 
* 

was chosen. 

Each MG process consists of three basic parts: relaxation, restric- 
tion, and interpolation [5]. 

The Relaxation Scheme 

The overall convergence rate of any MG process is greatly influenc- 
ed by the smoothing properties of the relaxation scheme. It can be shown 
analytically and experimentally that the usual multiple sweep marching 
[1] does not have good convergence and smoothing properties because short 
wave errors are not efficiently smoothed. Israeli and Lin [3] showed 
that certain modifications in the streamwise momentum equation, which 
vanish upon convergence, give rise to an iterative scheme which is equi- 
valent, in the linear case, to the SLOR method for one Poisson equation. 

In the general nonlinear case the modified iterative process is essential- 
ly equivalent to the relaxation of a single nonlinear Poisson-like equa- 
tion for the pressure. The velocities can be viewed as auxiliary vari- 
ables needed during the marching since they have no "memory" by them- 
selves. 

Furthermore, we have automatically gained the good smoothing 
properties of the line relaxation scheme of a single Poisson equation. 

The problems associated with the loss of ellipticity of the difference 

* Some of the elements of the present approach were used independently by 
Rubin and Reddy [7]. Detailed comparisons cannot be made because converg- 
ence rates and storage estimates were not presented there. 
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approximation for the Navier-Stokes equations at high Reynolds number 
[5] are thus avoided and no upstream-weighting or artificial viscosity 
are required. There results a considerable saving in storage, as well 
as a simpler relaxation scheme (compare to the distributive relaxation 
[5]) where the convergence rate is essentially independent of the 
Reynolds number. We note that the same marching algorithm can thus be 
used for the (subsonic) Euler equation with the same favorable converg- 
ence rate. (For supersonic flows the marching method is non iterative.) 


A part of the analysis of [3] is repeated here to motivate the 
later extensions to three-dimensional and compressible flows. We start 
with the PNS equations (2.1)- (2.3) and linearize them about a constant 
state. We also introduce 


L(f) = ■— (Uf) + (Vf) 
3x 3y 


_1_ 

Re 



f 


(4.1) 


where U and V are constant reference velocities. The next step is 
to discretize the equations only in the x direction to obtain: 


m-1 m 
Ax 


(V ) 
y m 


= 0 


D(U ) 
m 



m 


(4.2) 

(4.3) 


D(V ) = 
m 


-(P ) 
y m 


(4.4) 


where U , V and P are functions of y. Here D(f ) is the semi 
m m m m 

discretized form of E(f) at the marching station m. The semi - 
discretized system should be discretized also in the y direction before 
solution is attempted, but since the specific form of this discretization 
is not important for the following argument, we postpone this step for 
the sake of transparency. 


k k 

The marching iterative procedure assumes that U (y) , V (y) are 

k-1 ^ ^ 

known as well as P for m = 2,3,4,... M, where k is the current 

m 

iteration index. Therefore, the marching scheme for m £ 2 is: 


k k 
U -U 
m— 1 m 

Ax 


+ (V ) = 0 

y m 


P k - 1 -P k 

D(u k ) = _ 

m 


(4.5) 


Ax 


(4.6) 
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D(V k ) = 
m 


-<p k > 

y m 


(4.7) 


We now apply D to Eq. (4.5) and differentiate Eq. (4.7) with 
respect to y. Elimination of the V terms between Eqs. (4.5) and 
(4.7) gives: 


D(U k - U k ) =-(P k ) Ax 
m-1 m yy m 


(4.8) 


Now by substitution of Eq. (4.6) into Eq. (4.8) we get: 


P m+1 " " P k 1 + P k + Ax 2 ( pk ) = 0 . 

m+1 m m n-1 yy m 


(4.9) 


It follows that the marching scheme for the primitive system (4.2)- (4.4) 
can be viewed as a line iterative scheme for the semi-discretized Laplace 
equation; indeed upon convergence Eq. (4.9) will become: 


P m+r 2 V P m-l 


+ (P ) = o . 

yy m 


(4.10) 


In order to find out the rate of convergence of Eq. (4.9) to the final 
state (4.10), we Fourier transform Eq. (4.9) in y assuming appropriate 
boundary conditions in that direction: 


k r I ny k-1 „ I ny 

P = Z e J , P = Z e 1 
mm mm 


(4.11) 


where I = -1 and n is the Fourier wave number. After substituting 
these definitions into Eq. (4.9) we get: 


m+1 


- Z - 
m 


m 


+ Z 


m-1 


a 2 2- 

- Ax n Z = 
m 


(4.12) 


Transforming in the x direction we define again 

„ * I 

Z = Ae (4.13) 

where 0 is the wave number in the x direction. By substitution of 
this definition into Eq. (4.12) we get: 


l-l 

1 A * 


1-e 


-10 


' A 2 2 —10 
1+Ax n -e 


(4.14) 


This means that all the long waves (with small 


Ax n ) in the cross flow 
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direction are only weakly damped irrespective of their structure in the 
marching direction. In particular the n=0 modes which exist for 
derivative boundary conditions in the cross flow direction are not 
affected at all by the relaxation, i.e.. 



On the other hand, the well known SLR scheme for the Poisson equation 
gives (after the same Fourier transformations) : 


Z , - (2+ Ax 2 n 2 ) Z + Z = o 
m+1 m m-1 


(4.15) 


and 


I -I 

'a 1 


q-e 


-ie 


2 2 

• q = 2 + Ax n > 2 


and also 

|Ai2 


1 A 1 2 

q +l-2qcos0 


(4.16) 


(4.17) 


This quantity is less than 1 for all acceptable q's and cos0 < 1. 

Most waves are strongly damped* and only the longest waves in both direc- 
tions are weakly damped by the iteration. This behavior was used to 
accelerate the convergence as is done by the SLOR technique, Chebychev 
acceleration, or Multi-Grid method. 

The question is how to generate an equivalent relaxation scheme 
for the primitive system in the marching form. This means that we may 
add terms which can be evaluated during the marching process but should 
vanish upon convergence. 

A rational approach to the construction of the relaxation scheme 
is to retrace backwards the steps of the derivation of the discrete 
Laplace equation from the discrete primitive equations. We start from 
the SLOR equation (4.15): 


m+1 


- 2 2 - 
2Z + Z „ = Ax n Z 
m m-1 m 


which we inverse Fourier transform with respect to y to get: 


k-1 _k k A 2 k . 

P , - 2P + P _ = -Ax (P ) 
m+1 m m-1 yy m 


*It was pointed out by J. South that (4.9) can be viewed as an over- 
relaxed version of (4.10) with an over-relaxation factor: ^ = 2 which 
is not a good choice for £7. 
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Now, we substitute Eq. (4.8) for the right hand side of the last equa- 
tion to get: 


E ^' 1 

m+1 


2P k + P k l 
m m-1 


AxD(U k - U k ) 
m-1 m 


which can be written as: 


(Pi' 1 - pi) + (p! _ 1 - P*>- (P^- 1 - P k ) =AxD (U - u ) ; 

m m m-1 m-l m 


m+1 


m 


m 


adding the equations from m = 2 and using the linear form of D we get: 
-(P^-PV f (P k " : 


i=2 


P?) + (P k 1 ~P k ) =AxD(U -UJ 
1 ^ 1 ml 


but from Eq. (4.3) : 


P k 1 - P k = -AxD(u k ); 


therefore, we get for m > 2 


(4.18) 


D(U ) = - 
m 


p k ~t pk 

m+1 m 


Ax 


- h X ‘-*' 1 - > k > 


Ax . _ 
1=2 


(4.19) 


Eq. (4.19) contains all the modifications required in order to 
convert the iteration scheme of (4.5) — (4.7) into a scheme equivalent to 
the SLOR scheme for one Laplace equation with "over-relaxation" factor 
ft - 1 . We see that in this approach only the x momentum equation is 
modified. The new added term can be generated easily during the 
marching process and is inexpensive in storage (one extra line vector) 
and computation (one substruction per grid point). In what follows we 
derive Eq. (4.19) in a more general way and introduce the over- 
relaxation parameter ft > 1. 

In practice we will use difference approximations and boundary 
conditions also in the y direction, and the resulting scheme may not be 
amenable to the discrete analogue of the Fourier transform. It is 
therefore worthwhile to generalize the previous approach by using the 
matrix finite difference formulation. 

Let the vectors U^, V , P^ contain the N values of the cor- 
responding variables on the m-th line (x = constant) of the marching 
sweep (including the specified boundary values). The U-momentum 
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equation (1.5) can be written in the form: 


p i i - P = -AxD (U ) 5 R . 
m+i m m m 


(4.20) 


On the other hand elimination of V between the continuity and the 

m 

V-momentum equations will result in: 


FP = R 
m m 


b m-1 


(4.21) 


2 3 

where F = Ax I — j . Substracting successively U-momentum equations 
(4.20) and using^ Eq. (4.21) gives: 


P - (21+ F)P +P = 0, m = 3,4,... 
m+i m m-1 


(4.22) 


which is Laplace's equation. The first equation of (4.20) can be used 
as a derivative condition at the left (inlet) boundary, namely: 


P 3 " P 2 " R 2 


(4.23) 


We now apply the SLOR scheme to the last two equations (ignoring 
temporarily the downstream boundary condition) to get the downstream 
marching form: 


< * - r 2 
Jk) 


' m-1 


(21 + F ) P* + P^“ 1} =0, m = 3,4,... 
m m+1 


(4.24) 

(4.25) 


(k) * (k-1) 

where P = ftp + (l-ft)P ; ft is the overrelaxation factor, and 

m m m 

the superscript denotes the iteration sweep number. In order to recover 
the primitive variable formulation, we relate the velocity field in Eq. 
(4.21) to the starred pressure field, i.e.. 


FP = R - R „ , 
m m m-1 


Substitution in Eq. (4.25) gives: 


P (k ^ - 2P + P (k 1) = R - R , m = 3,4,, 
m-1 m m+1 m m-1 


(4.27) 


Successive summations of Eqs . (4.24) and (4.25) give: 

(k-1) * 

- p = * + ' m = 2,3,4,... 

m+ 1 m m m 


(4.28) 
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which is the primitive variable marching form of the U-momentum equation. 


The source term S m in Eq. (4.28) satisfies: 


S = S . + (p 
m m-1 m 




(4.29) 


with 0. It can be seen that S^ vanishes upon convergence The 

computational form of (4.28) for m = 3,4,... is: 


* 


-2P 

m 


R 

m 


+ S 

m 


(4.30) 


S = 
m 



(k-1) 

m+1 


+ 2P 


m-1 


(k) 

m-1 ' 


* 




(4.31) 


Thus, the theory of overrelaxation can be applied exactly to the 
constant coefficient case of system (2 . 8) - (2 . 10) . For the non-linear 
case this theory can serve as a guide to the choice of Q. Alternately, 
one can choose H = 1 and apply the Multi-Grid procedure. 


Restriction and Storage Requirements 

Let the finite difference approximation of equations (2.1)- (2.3) 
on the finest grid M be represented as in [5] : 


iWfx) = Fj (x) 


(4.32) 


. MMMT 

where x - (x,y) , W = [U ,V ,P ] is the exact solution of the dif- 
ference equations, and j is the number of the differential equation, 
j = 1,2,3. 

The problem is transferred from the current level k to a coarser 
level k-1, see Fig. 4, by correcting the right hand side of (4.32) 


k-1 - T k-1 ,~k-l~k . 

F j 1x1 ' b (I j,k v ' !x!! 


, T k-1. k,~ _k~k,~., 

+ v If . (x) - l w (x) ] 
J J J 


(4.33) 


in the Full Approximation Storage (FAS) mode. w k (x) is an approximation 

to W (x) in the finer level. X . ^ and i k ^ are proper restriction 

D J t k 

operators for equation j. 

The term in square bracket in equation (4.33) is the residual of the 
j-th equation. For the present marching scheme there is no residual in 
the continuity and in the y-momentum equations since they are solved 
exactly in each step. The residual of the x-momentum equation results 
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only from the streamwise pressure gradient term and its computation needs 

~k-l 

only one substraction . I. , 

3i k -i ^k-i~k - k-i 

which yields for the continuity equation: L (I_ w (x) ) = 0. I 

1 1 ,k j ,k 


was chosen to be linear interpolation, 


j = 1/2 is computed by averaging in both the x and y directions. 

k— 1 ...... 

1^ ^ is a simple injection. 


In summary, equation (4.33) takes the following terms: 
k-1 . ~ . 

F (x) = 0 

k-1 ~ _ k-1 -k-l~k - k-1 k r k-k ~ 

F 2 ( 5 L 2 (1 2,k W I 2,k (F 2~ L 2 W (x) } 

_k-l - k-1 -k-l-k ~ 

F 3 (X) = L 3 (I 3,k W (X)) ' 


Two consequences should be emphasized: 

k-i - k-1 - 

(a) Only two corrections (F^ (x) , F 3 (x) ) have to be computed and 

stored. 

(b) All the dependent variables must be transferred in order to compute 

k-1 ~k-l~k ~ 

the corrections (L. (I. w (x) ) , j = 2,3). Since only the pressure 

J J fk 

is stored, these corrections must be computed during the marching 
process . 


It follows that in addition to the pressure on all grids, one has 
to save one correction term for each momentum equation on the coarser 
grids. Assuming N computational points on the finest grid, a simple- 
minded estimate gives 3SN/7 storage locations for the three-dimensional 
NS Multi-Grid solution, and 11N/7 for the PNS marching MG solution. 

For the two-dimensional case the corresponding figures are 14N/3 and 
6N/3. 


Interpolation 

Since the present marching scheme generates the velocity field from 
the pressure, only the correction to the pressure must be interpolated 
back to the fine grid. 
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5. GENERALIZATIONS 


In order to generalize the preceding approach we note that the 
essence of the relaxation procedure is the replacement of the term 
Ax3P/9x by the marching difference form: 


4x|£ 

dX 


-(2P* + S ) 
m m 


(S.l) 


where is already known. If (5.1) is differenced it will (using 

the definition of S ) give rise to 
m 3 


2 (P 


P ,) + S 
m-1 m 


— V- 1 * v 

s . = -(P* 7 - 2P + P* ) 
m-1 m+1 m m-1 


Thus, the correct successive line over relaxation form is implicitly 
obtained for the second derivative of the pressure in the marching 
direction. (It should be emphasized again that the second order elliptic 
equation for the pressure is neither derived nor used in the algorithm 
itself. ) 

The implication of the present technique is: if it is known that 
the equations can be manipulated so that some variable will satisfy 
approximately a second order elliptic equation, we should use the replace- 
ment (5.1) for the derivative of that variable in the main flow direction. 
An efficient marching scheme will thus be generated. 

The present version of the algorithm will be applied to the sub- 
sonic compressible multi-dimensional Navier-Stoke ' s equations. Several 
particular cases will be examined. 

The first step is the derivation of an elliptic equation starting 

with : 

V*VV=-i Vp + vV 2 V + V(V*V) (5.2) 

where V is the velocity vector. In addition we will require the equa- 
tion of state of a perfect gas 

p = PRT (5.3) 

and the continuity equation in the form 

V*V = -V*V£np. (5.4) 


It follows from (5.3) that 
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i Vp = RTV£np + RVT 
and therefore 

Vp = RTV 2 £np + RV 2 T + l.o.t. 

(l.o.t. stands for lower order terms) . 

Also 

V-(V-VV) = (V-V) (V-V) + l.o.t. = ~(V-V) 2 £np + l.o.t. 


(5.5) 


(5.6) 


(5.7) 


Taking the divergence of (5.2) and using (5.6) and (5.7) we get 

(V*V) 2 £ n p = RTV 2 £np + RV 2 T + y V 2 (V-V) + l.o.t. (5.8) 

Several special cases follow: 

1) Incompressible case . Here V*V = 0, and we get directly (taking the 

divergence of (5.2)): 

v 2 p = l.o.t.; (5-9) 

thus the gradient term was differentiated once and the replacement 
(5.1) should apply in two or three dimensions . In the later case 
we have to compute simultaneously all the variables in the marching 
plane., m, and so we get a successive plane overrelaxation scheme. 
It is possible that an alternating direction scheme can be used to 
solve the coupled system in the m th plane, but a multi-grid 
approach seems to be preferable. At the present time numerical 
results for the three dimensional case are not available. 

2) Isothermal case . We get 

2 

(V*V) 2 £np = V 2 £np + l.o.t. + v.t. 

where a 2 = yRT is the adiabatic speed of sound, y is the ratio 
of specific heats, and v.t. is the viscous term to be discussed 
later. 

3) Isentropic case . Here, (y-l)V£np = V£nT and therefore: 

(V*V) 2 £np = a 2 V 2 £np + l.o.t. + v.t. 

4) Constant stagnation enthalpy . Here, a 2 = a 2 + • v 2 and we get 

(V* V) 2 £np + (y-l)V* (V*V) -V = a 2 V 2 £np + l.o.t. + yv.t. 

In all the compressible cases considered the prominent balance is: 
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(V* V) 2 £.np = a 2 V 2 <»np . 

Using local stream aligned coordinates s and n, we find 
(V-V) 2 = V 2 

8s 2 

V 2 = _3i + ; 

3s 2 9n 2 


therefore p = £np appears in the form: 


(1-M 2 ) -■ — + — - = other terms. 
9s 2 9n 2 


After the parabolization of the viscous term, only the left hand side 
has a second derivative of p in the streamwise direction. 
Specifically W 2 V*V is replaced by 


9n 2 



V*V£np . 


This term cannot become large since the pressure does not have large 
gradients in the boundary layer. 

We argue that if our iteration is appropriate for the p equation 
it will be a good scheme overall. 

To get a successive line (or plane) over relaxation scheme, all we 
have to do is replace all the occurrences of ^ with the marching form 
(5.1). All the properties of Section 4 will be the same as long as 
M 2 < 1. 

In fact, better convergence can be expected as M 2 approaches 1 
since the quantity q of (4.16) will become now 2+ (Ax 2 n 2 /1-M 2 ) . 

Only p will have memory and must be globally saved and updated by the 
iteration procedure, p will also transmit the downstream information 
and must be specified there. 

For transonic flows a conservation form is preferred and it may be 
more convenient to work with p rather than £np . An elliptic equation 
can be derived for p, but care must be taken to transmit the downstream 
information via p. Upstream information should not be transmitted by 
p and p should not be specified at the inflow; otherwise, the problem 
will be overspecified. Consider, for example, the term 9pu 2 /9x ; it is 


discretized as 
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P ,U2 
m-1 m-1 


) 


however, at the station m we should compute the U m velocities coupled 
with the P m _-L densities. The approach of Reddy and Rubin [12] where the 
pressure is specified both at inflow and at outflow is inconsistent unless 
one happens to know the right pressures before the computation. The 
inconsistency and consequent error can be easily demonstrated by one- 
dimensional examples. 


6. RESULTS 

In order to check the MG algorithm, we choose the following analytic- 
al solution. It satisfies the continuity equation but gives rise to 
source terms in the momentum equations: 

U = A + (x+y)™; V = -(x+y)™; P = -(E1+E2) (x+y)™ (6.1) 

where a and b from equations (3) are defined by: 

a = El + F(x+y) n ; b = E2 - F(x+y) n (6.2) 

and El = 1? E2 = .2; F = .2; A = 5; Re = 1000; m = 4; n = 2. The 
coarsest grid consists of 4x4 intervals. 

Figure 5 compares the MG convergence history of different relaxation 
schemes. In the MG solutions three levels were involved (M=3) . The 
horizontal coordinate gives the number of Work Units (WU) , where each 
work unit is equivalent to one global iteration on the finest grid. The 
vertical coordinate gives the logarithm of the dynamic residual e. The 
dots show the solution of the equivalent Poisson equation (with the same 
solution for the pressure but with Dirichlet condition over all the 
boundaries) . The linearized PNS equations were solved with and without 
the streamwise pressure gradient correction of [3]. The corresponding 
(17 x 17 points) single grid convergence history is plotted for compar- 
ison (for the case of Q = 1). The corrected discrete equations and 
the Poisson equation exhibit very similar convergence whereas the conver- 
gence of the unmodified equations is much worse. Upon increasing the 
number of grids in the unmodified equations, the convergence deteriorates. 

The Reynolds number independence of the scheme is demonstrated in 
Figure 6, where the convergence history is presented for Reynolds numbers 
1, 10 3 and infinity. 

In order to check the non-linear version of the code, several test 
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cases were run; the incompressible flow over a flat plate, the flow 
along an axisymmetric cylinder, entrance flow between two flat plates, 
and the flow behind the trailing edge of a flat plate. In all cases 
good agreement was obtained with known solutions. The details will be 
presented elsewhere. Here we show (Figure 7) the convergence history 
for a flow over a flat plate with uniform upstream profile and Neumann 
condition for the pressure at the exit. While the number of levels is 
varied, the finest grid remains the same and consists of 65 x 65 points. 
In Figure 8 there is a comparison between the present results for the 
flow near the trailing edge of a flat plate and the results of refer- 
ence [11]. The skin friction coefficient CF is shown for z < 1 while 
the center line velocity UC is shown for z > 1. The trailing edge is 
at z = 1. 
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Figure 1: Pressure Residue vs. the Global Iteration Number. 


B 

UPSTREAM 
any two of (u,v,p) 

A 


any two of(u,v,p) 

^ » 0 

DOWNSTREAM 

only (p) 

D 

any two of(u,v,p) 


Figure 2: Example of Permissible Boundary Conditions. 
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Figure 4: Relative placement of variables on two successive grids 
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